Phases of the one-dimensional Bose-Hubbard model 



T. D. Kiihner and H. Monicn 
Physikalisches Institut der Universitdt Bonn, D-53115 Bonn, Germany 
(February 1, 2008) 

The zero-temperature phase diagram of the one-dimensional Bose-Hubbard model with nearest- 
neighbor interaction is investigated using the Density-Matrix Renormalization Group. Recently 
normal phases without long-range order have been conjectured between the charge density wave 
phase and the superfluid phase in one-dimensional bosonic systems without disorder. Our calcula- 
tions demonstrate that there is no intermediate phase in the one- dimensional Bose-Hubbard model 
but a simultaneous vanishing of crystalline order and appearance of superfluid order. The complete 
phase diagrams with and without nearest-neighbor interaction are obtained. Both phase diagrams 
show reentrance from the superfluid phase to the insulator phase. 
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Quantum phase transitions in strongly correlated sys- 
tems have attracted a lot of interest in recent years. Usu- 
ally the basic particles are electrons, but in some inter- 
esting cases the relevant particles are not fermions but 
bosons. Examples of experimental systems with super- 
fluid and insulating phases^are Cooper pairs in thin gran- 
ular superconducting filmsEJ and cooper pairs or fluxes in 
Josephson junction arraysB. While injdimensions greater 
than one the existence of supersolidsO, i.e. phases with 
simultaneous superfluid and crystalline order, has been 
established in theoretical work, the situation in one di- 
mension is less clear. Recently normal phases that are 
neither crystalline nor superfluid have been found m a 
one-dimensional model of Josephson junction arraysQ in 
the region where supersolids are found in higher dimen- 
sions. In this paper we will verify whether supersolids or 
normal phases exist in the more general Bose-Hubbard 
model in one dimension. 

The Bose-Hubbard model contains the basic physics of 
interacting bosons on a lattice. It is a minimal bosonic 
many-particle model that cannot be reduced to a single 
particle model. The bosons have repulsive interactions, 
and they can gain energy by hopping to neighboring sites 
on the lattice. The Hamiltonian with on-site and nearest- 
neighbor interactions is 

HbH - -tY.^{b\h,+l + br^+l) - 

+ [/ n,{n, ~ l)/2 + VY.^ ri,n,+i , (1) 

where hi are the annihilation operators of bosons on site 
i, hi — b\b^ the number of particles on site i, t is the hop- 
ping matrix element. U and V are on-site and nearest- 
neighbor repulsion, and ^ is the chemical potential. The 
energy scale is set by choosing U = 1. 

The range of the interactions depends on the individual 
experimental situation. In general the lattice underlying 
the system is not an atomic lattice, but a larger structure 
like a Josephson-junction or a grain in a superconductor. 
In Josephson-junctions the relevant bosons can be cooper 
pairs or fluxes, resulting in different interactions. 

As a starting point we first consider the case of on- 
site repulsion only. In the (/i, t)-plane Mott-insulating 



regions are surrounded by the superfluid phaseB. These 
phases are separated by two types of phase transitions. 
On the constant density line the transition is driven by 
phase fluctuations and is of the Berezinskii-Kosterlitz- 
Thouless (BKT) type. The phase transition at the sides 
of the insulator, the generic phase transitiodj, is driven 
by density fluctuations. 

The Mott-insulating phases have integer densities and 
are incompressible, at the generic phase transition to the 
compressible superfluid phase the density of the system 
changes from commensurate to incommensurate. The 
characteristic energy Eg^^'^ of this transition is the en- 
ergy it costs to create a particle (p) or hole (h) exitation 
in the-, system. To calculate this energy we use defect 
statesO with the density of the Mott-insulator plus one 
additional particle or hole. Since the defect states and 
the insulator groundstate have fixed densities, a change in 
the chemical potential by A/i does not change the states 
themselves, but shifts their energy by A/^iN, where N 
is the total particle number. Taking into account that 
the characteristic energy is zero at the phase transition, 
this gives Ef''\^x,t) =\ fj.itfc'^''^ - ^ , where ^^(i) is 
the chemical potential at the phase transition and the 
critical exponent zv = IB. 

We use the infinite-size algorithm of the density-matrix 
renormalization group (DMRG)lI with periodic bound- 
ary conditions to determine Eg^'^^ . While the maximum 
number of particles per site in the Bose-Hubbard model is 
n = cx), it has to be cut off for practical calculations with, 
the DMRG. A maximum occupation number of n = 5u 
turned out to be sufhcicnt. 

Since the two defect states, one with an additional par- 
ticle, one with an additional hole, are needed to calculate 
the energy Eg''^^ , they are used as additional target states 
in the DMRG. Systems of up to 76 sites are calculated, 
keeping 128 states in each iteration. The chemical poten- 
tial fic'"'^\t) —\ El^^\t) — ^{t) I of the phase transition 
is calculated for various system sizes. The thermody- 
namic limit is found by extrapolating to infinite system 
size (Fig.|l|) . Repeating this calculation for various t gives 
the phase boundaries. 
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FIG. 1. The particle and hole excitation gaps plotted 
against the inverse system size. The straight lines are fitted 
to find /ic in the thermodynamic limit. The same scaling 
behavior has been found in all cases, {p = 1, t — 0.25, U = 1, 
V = 0.4) 



At the tips of the insulator lobes the phase transition 
is driven by phase fluctuations. The characteristic length 
of this transition is the correlation length 



r r 

of the correlation function 

r(r) = (6t(0)6(r)). 



(2) 



(3) 



The corresponding energy Eg is the energy gap between 
the groundstate and the first excited state with the same 
density: Eg ~ with the critical dynamical exponent 
z = l. 

The correlation length in the thermodynamic limit can 
be found by extrapolating from finite systems: — 
Coo + cl/L + b exp(— L/c), where L is the system size 
and a, b, c are fitting constants. The exponential term is 
small {a/b > 2), with c w 3, the results are not changed 
significantly by neglecting the exponential term. The 
phase transition is of the BKT type, where the corre- 
lation length diverges like 



const. 
exp( ^ =) 



t' 



(4) 



One way to findj-the critical point is fitting ^ to the 
calculated dataflilil. But by changing the fitting param- 
eters, this function can go to zero arbitrarily slow, hence 
this method is very sensible to numerical errors and the 
choice of data points, and we will not use it. 

Instead we locate the BKT transition using thci-anal- 
ogy of the superfluid phase to the Luttingcr liquidll3. In 
the superfiuid phase the correlation function T decays 
algebraically: 



r^-^(r) cx r 
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FIG. 2. K plotted against t. Kc = 1/2 at tc = 0.277 ±0.1 
([/=!, l/ = 0). 

The exponents Kc at the ukase transitions are known 
from Luttinger liquid theoryll3. At the BKT transition 
with p = 1, if is expected to be Kc = 1/2. An alge- 
braic function is fitted to the correlation functions cal- 
culated with DMRG for different system sizes. Due to 
the periodic boundary conditions the decay of T is very 
close to algebraic even in small systems. The thermo- 
dynamic limit of K is found by extrapolating K{L) — 
K — const./ L. Since the decay of r'^-'^(r) is very close to 
(^, the main source of errors in K is this extrapolation 
from finite system sizes. We find the phase transition at 
tc = 0.277 ±0.01. 

This is in good agreement with tc — 0.2X5 ± 0.005 
found in an exact diagonalization approaclxl, and in 
qualitative agreement with the Bethe-Ansatz solution 



BA 



1/(2^3) « 0.289 for the trut 



lated model with a 
Three works find. 
tc = 0.304±0.002EI 



maximum of n = 2 particles per s: 

somewhat bigger valueSpL — 0.298I , ^ _ _ _ _ _ _ 

and t = 0.300 ± 0.00313. r-Parly QMC simulations re- 
sulted in tc = 0.215 ± 0.0 iN. The range of these results 
demonstrates that determining the location of the BKT 
transition is ill-conditioned. 

Fig. ^ shows the phasediagram in the (/x, t)-plane, in- 
cluding the generic phase boundaries and the location of 
the BKT transition. For t > 0.225 we find that the lower 
phase boundary is bending down. This means that the 
Mott-insulator phase is reentrant as a function of t, an 
unusual feature that has not been observed before. It im- 
plies that increasing t, which corresponds to increasing 
the kinetic energy, can lead to a reentrance phase transi- 
tion from the superfluid phase to the insulator phase. 

In a study inspiijed by this work a high order strong 
coupling expansionim was used to determine the phase 
diagram. The phase boundaries found in that work are 
in excellent agreement with our DMRG results, demon- 
strating the high numerical accuracy that can be achieved 
with the DMRG. The strong coupling expansion study 
confirms the existence of the reentrant phase transition 
first found in this work. 
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FIG. 3. The phase diagram without nearest-neighbor in- 
teraction (MLMott-insulator with density one, SF:superfluid 
phase). The solid hnes showi-ai Fade analysis of 12th or- 
der strong coupliagi expansiondlZI, the boxes show Quantum 
Monte Carlo datsO. The circles are the DMRG results, the 
dashed lines indicate the area with integer density. The er- 
ror bars in the ^ direction are smaller than the circles, the 
error bar in the t direction is the error of the BKT transition. 
(f/ = 1, y = 0). 

One experimental realization of a one dimensional 
bosonic lattice system is provided by fluxes, in a 
Josephson-junction array. In a recent experimentB Mott- 
insulators with flux densities 1/3, 1/2, 2/3, 1.. were 
found. In the Bose-Hubbard model insulators with these 
densities can only appear in the presence of longer ranged 
interactions. A first step in understanding these longer 
ranged interactions is studying the effect of repulsion be- 
tween nearest neighbors. 

In the presence of nearest-neighbor interaction a new 
insulator phase appears at half integer densities. It is a 
charge density wave phase (CDW) with a wavelength of 
two sites, and like the Mott-insulator at integer density 
it has an excitation gap and is incompressible. The crys- 
talline order is characterized by the non-zero structure 
factor 

in the thermodynamic limit is determined by extrap- 
olating from DMRG calculations for finite systems. A 
maximum particle number of four particles per site is 
chosen for p — ^. Since the groundstate and the first 
exited state in the CDW are degenerate in the thermo- 
dynamic limit, and close to degeneracy in finite systems, 
the first exited state is used as an additional target state. 

is found to scale like 5'^(L) = STr+a/L+b exp {—c/L). 
The exponential term gives only a very small contribu- 
tion {a/b > 10), c is of the order 10 - 20. 

In contrast to the transition from the Mott insulator 
to the superfluid, which is governed by superfluid order 



only, the transition from the CDW to the superfluid is 
governed by superfluid and crystalline order. 

There are three possible scenarios for this transition: 
a) There is a direct phase transition - the vanishing of 
crystalline order and the appearance of superfluid order 
coincide, b) There is an intermediate phase with simul- 
taneous superfluid and crystalline order, the so-called 
supersolid phased, c) There is an intermediate normal 
phase with neither superfluid nor crystalline order. 

In higher dimensional bosonic systems supersolids ex- 
ist, but they ]myc not been observed in one-dimensional 
systems so fartl. Recently a normal phase (scenario c) 
was found in a numerical study of the one-dimensional 
Quantum-Phase modelQ, which is the high density limit 
of the Bose-Hubbard model. This raises the question 
whether such a normal phase also exists in the Bose- 
Hubbard model. 

The correlation length ^ characterizing superfluid 
order diverges in the superfluid phase. If the structure 
factor S-jr is also governed by this correlation length, there 
is a direct phase transition from the CDW phase to the 
superfluid phase. In that case, close to the phase tran- 
sition S„ ~ C^^^i^/L), where $ is a scaling function. 
Note that this functional form cannot be transformed to 
a power law behavior depending on t due to the BKT be- 
havior of ^ (Q) . To verify the existence of <I> at the CDW 
to superfluid transition, we plot S*^ versus ^ (Fig.^). We 
find Stt ~ ^-0-4±o.i^ which shows that there is a direct 
phase transition and no normal or supersolid phase. 
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FIG. 4. Stt plotted against ^ at the BKT transition of the 
charge density wave phase {U = 1,1^ = 0. 4). The slope is 
0.4 ±0.1. 



The superfluid stiffness in the Luttinger Liquid is al- 
ways non-zero, even for large KB. It ha s be en shown 
that for K > 1 a single weak link or barrieill3lH3 becomes 
relevant, reducing the superfluid stiffness to zero by effec- 
tively cutting the system in two parts. ,Since K = AdX the 
sides and X = 2 at the tip of the CDWll3, this means that 
the CDW is surrounded by a region with K > 1. The 
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normal phase found by Baltin and Wagenblasto was ob- 
served in a Quantum Monte Carlo study at finite temper- 
atures, where the zero-temperature phase diagram was 
extracted with finite-size scaling. Baltin and Wagenblast 
suggested that the normal phase might be the region of 
the Luttinger Liquid with A' > 1. But since their cal- 
culations were for a system without impurities, which at 
incommensurate densities always has a superfluid stiff- 
ness, they should not have observed an effect caused by 
an impurity. A possible explanation of their result is that 
thermal fluctuations might have a similar effect as a weak 
link. This is supported by the fact that Baltin and Wa- 
genblast could not determine the scaling function for the 
dependence of the superfluid stiffness on system size and 
temperature. 



tentials, which shows the reentrant behavior already ob- 
served in the case without nearest-neighbor interaction. 

In conclusion, we have presented methods to determine 
the generic as well as the BKT phase transitions of the 
Bosc-Hubbard model with the DMRG. At the tips of the 
insulating regions we found a reentrant phase transition 
from the superfluid phase to the insulator. Including 
nearest-neighbor interactions we obtained the new phase 
diagram and demonstrated that there is no normal or 
supersolid phase, but a direct phase transition from the 
CDW to the superfluid phase. 

We would like to thank T. Giamarchi, A. J. Mil- 
lis, R. Noack, A. v. Otterlo, G. Schon, H. Schulz and 
S. R. White for useful and interesting discussions. 
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FIG. 5. The phase diagram of the Bose-Hubbard model 
with nearest-neighbor interaction (MI: Mott-insulator with 
density one, CDW:charge density wave with density one half, 
SF:superfluid phase). The error bars in the ^ direction are 
smaller than the circles, the error bars in the t direction indi- 
cate the errors of the BKT transitions. {U = 1, V — 0.4). 



The onset of superfluidity is again determined by the 
decay of the correlation functions. At the CDW the criti- 
cal exponent is K^'^^ = 2. The BKT transition is found 
at {t/U)^^^ = 0.118±(m04. This is in agreement with 
tc « 0.1 found with QMCB. For the Mott-insulator with 
density p = 1 the exponent K is found to change very 
slowly close to the phase transition, causing a high er- 
ror margin in our calculation. We find the critical value 
Kc = 1/2 at tc w 0.325 ± 0.05. This indicates that the 
critical point is shifted to higher ratios of t by increasing 
V. Within the numerical accuracy the critical point may 
also be independent of V. This cor^adicts QMC results 
that tc is reduced if V is increasedtj. For V = 0.4 they 
found tc ~ 0.17. 

Fig. 1^ shows the phase diagram of the Bose-Hubbard 
model with nearest-neighbor interaction in one dimen- 
sion. To our knowledge this is the first time this phase 
diagram has been calculated. The tips of the insulating 
regions are bending down towards smaller chemical po- 
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